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We study the ground-state properties of the two- 
dimensional quantum spin systems having the valence-bond- 
solid (VBS) type ground states. The "product-of-tensors" 
form of the ground-state wavefunction of the system is uti- 
lized to associate it with an equivalent classical lattice sta- 
tistical model which can be treated by the transfer-matrix 
method. For diagonalization of the transfer matrix, we em- 
ploy the product-wavefunction renormalization group method 
which is a variant of the density-matrix renormalization group 
method. We obtain the correlation length and the sublattice 
magnetization accurately. For the anisotropically "deformed" 
5 = 3/2 VBS model on the honeycomb lattice, we find that 
the correlation length as a function of the deformation pa- 
rameter behaves very much alike as that in the S — 3/2 VBS 
chain. 



I. INTRODUCTION 

There are many applications of the density matrix 
renormalization group (DMRG) (l]j|] which was origi- 
nally applied to the one-dimensional (ID) quantum sys- 
tems Jl,§. Due to its remarkable success, the DMRG has 
now become one of the standard methods for studying ID 
quantum models, two-dimensional (2D) classical mod- 
els H and (1 + l)-dimensional classical non-equilibrium 
models §§. 

As was pointed out in Refs. ||[7J, DMRG is a vari- 
ational method under the matrix-product-form ansatz 
(MPFA) for trial wavefunctions whose usage dates 
back to the work of Kramers and Wannier [p|dl0|. 
This point of view leads to some non-trivial refor- 
mulations of the DMRG: the direct variational ap- 
proach, |6p|, pd| , [l2]| the product-wavefunction renormal- 
ization group (PWFRG), [jl3| the corner-transfer-matrix 
renormalization group (CTMRG). §JlO) Further, in this 
view, the success of the DMRG implies the unexpected 
accuracy of the MPFA wavefunctions. 

Having seen the success of the DMRG for ID quantum 
and 2D classical systems, it is natural and important to 
explore its higher-dimensional generalizations. For 2D 
quantum case, the ladder approach, [|l4|-|l6|] can be re- 
garded as such, but it is essentially the one-dimensional 
algorithm. For 3D classical case, as an extension of 
the CTMRG, |9|JlCt| the "corner-tensor" approach has re- 
cently been proposed, Jl7[ which is the first among the 
truly higher-dimensional algorithms, although it may not 
be the only possibility. 



In our view, the most promising one in generalizing the 
DMRG to higher dimensions seems to be the one which 
is based on a generalization of the MPFA. We should 
remark that, we have already known a one: a wavefunc- 
tion which is a product of local tensors, typical example 
being the valence-bond-solid (VBS) state ]l8| , |l9t for a 
class of 2D quantum antiferromagnets. It has been well 
known that the wavefunction of the ID VBS state can 
be expressed as a product of 2 x 2 matrices, implying 
that the MPFA is exact for the system; it has also been 
known that the wavefunction of the 2D (or higher dimen- 
sional) VBS state can be expressed as a product of "ten- 
sors" (generalized objects of matrices). Let us simply call 
such form of wavefunction a tensors-product-form (TPF) 
wavefunction. A natural generalization of the MPFA is, 
then, the tensors-product-form ansatz (TPFA). Accord- 
ingly, we can formulate the TPFA-variational method 
which can be thought of a "higher-dimensional DMRG" . 
The aim of the present article is to make the first step 
in this TPFA-variational approach to higher-dimensional 
quantum systems. We focus on the property of the VBS- 
type state in 2D, and for this purpose, we develop a re- 
liable numerical method at the same time. 

The models we consider in the present article is the 
2D quantum antiferromagnets which have the VBS-type 
ground states; jl8-20| the TPFA is exact for these mod- 
els. Let |VBS) be the unnormalized ground-state vector 
whose wavefunction is exactly given by a product of local 
tensors. A main part of our problem is to evaluate the 
expectation value of a given observable O 



(O) = (VBS|C|VBS)/(VBS|VBS). 



(1.1) 



For ID VBS-type models, due to the matrix-product- 
form structure of |VBS), the RHS of (1.1) can be in- 
terpreted [|ll|fL9| as a thermal average in a ID classical 
statistical-mechanical model; the transfer-matrix method 



allows us to evaluate the RHS of (1.1) exactly. For 2D 
VBS-type models, similar interpretation as a 2D classi- 
cal statistical-mechanical problem is straightforward due 
to the TPF structure of |VBS). Based on this inter- 
pretation, Niggemann et al. made the classical Monte- 
Carlo study of the anisotropically generalized S = 3/2 
VBS model on the honeycomb lattice. The approach 
we take in the present work is, in a sense, a direct gen- 
eralization of the ID case; we treat the associated 2D 
classical statistical-mechanical problem by the transfer- 
matrix method. What is essential in our approach is 
that, for diagonalization of the transfer matrix, we em- 
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ploy the DMRG (to be precise, PWFRG, in the present 
article) allowi ng u s to make highly reliable, close-to-exact 
evaluation of (1.1). The use of the DMRG also has an im- 
portant implication, in the light of the TPFA- variational 
method: 2D TPFA-variational calculation is reduced to 
ID TPFA-variational method, namely the DMRG. Ac- 
cordingly, we can, in principle, formulate the "nested" 
TPFA-variational approach where Z?-dimensional TPFA- 
variational calculation is reduced to (D — l)-dimensional 
one, which in turn is reduced to {D — 2)-dimensional one, 
and so on. 

This paper is organized as follows. In Section |l[ we 
give explicit form of local tensors for the 2D VBS-type 
models. For actual calculations, we consider the 5 = 2 
isotropic VBS model on the square lattice, and the "de- 
formed" S — 3/2 VBS-type model on the honeycomb 
lattice proposed by Niggemann et al. In Section [II, we 



explain our calculation method. We extend the PWFRG 
so that it can handle the asymmetric transfer matrices 
which we encounter in treating the 2D VBS-type mod- 
els. In Section |vL calculated results are given. The last 
section (Section |Vj is devoted to summary and conclu- 



II. VERTEX MODEL ASSOCIATED WITH THE 
VALENCE-BOND-SOLID STATES 

A. Two-Dimensional Valence-Bond-Solid 
Wavefunction and Local Tensors 

Construction of the higher dimensional VBS state has 
already been made in the original paper. |l£| We shall 
give a brief explanation of the TPF structure of the gen- 
eral VBS state. 

Consider a spin-5 operator situated at a site i on a 
lattice which may not necessarily be a "regular" one. 
Let {|cr)} (ct = —S, —5 + 1, . . . , S) be the correspond- 
ing s z -diagonal basis set. We assume that the site i is 
"M-valent" with M = 25. We then prepare M spin- 
1/2 operators, and regard the spin-5 operator as the one 
constructed from these component operators. Let \r) a ) 
(rj a = ±1/2, a = 1, . . . , M) be the s z -diagonal base of 
the a-th component spin, and denote 

IM> = vm) = \m) ® \m) \vm)- (2.1) 

Consider the symmetrized state 

KM)) = I(?7i,?72,---,?7m)) = J\f n Sym n \r)i,r) 2 ,...,riM), 

(2.2) 

where Sym^ stands for the symmetrization with respect 
to the indices {rj a }, and Af n is a factor introduced for 
the proper normalization, (({?7})|({^})) = 1- Clearly, we 
have 



We then define a 
M<*\{v})(= A(a\r] U r] 2 ,...,T] M )) by 



local 



tensor 



A(v\{r,}) = (v\({r,})), (2.4) 
which is the building block of the VBS state (see Fig. [l]). 



fig-1 



Note that, by definition, the tensor A(a\rji, 772, . . . , tjm) 
is symmetric with respect to the indices {rja}. Further, 
we adopt the phases of the vectors {\a)} so that the ten- 
sor has the property 



A(a\{r,}) = A(-a\{-r,}). 



(2.5) 



Let Ai(ai\{ri 1 }) be the local tensor associated with the 
site i. The VBS wavefunction ^vbs(cI) C2, • • • , vn) is 
given by 



i) <hj> 



(2.6) 



where we have introduced the sign factor $y (rf , rj 3 ) for 
each connected site pair < i,j >, corresponding to forma- 
tion of the "valence bond" . [ |18||l9| For the valence bond 
formed from the a-th bond at the site i and /3-th bond 
at the site j, the sign factor is explicitly given by 



where 



S ij ( V \rf) = e( V i,4), (2.7) 



1 (f = -1/2 and £' = 1/2) 

-1 (1 = 1/2 and £' = -1/2) . (2.8) 

(otherwise) 



The sign factor can be absorbed into the definition of 
modified local tensor as 

A(a\r)x, ...,fj a ,...,) = ^2e{f] a ,r]' a )A(a\T)i, ...,»£,...,) 

(2.9) 

We can introduce modified tensors where two or more in- 
dices (at a-th, /3-th, . . ., positions) are replaced by those 
with bars (A(a\r)i, . . . ,fj a , . . . ,fjp, . . .), etc.). In terms of 
the modified tensors, we can express the VBS wavefunc- 
tion as a product of these tensors summed over the bond- 
variables {77} in a form similar to (2.6), with the sign fac- 
tor Sij(ri % ,rf) simply replaced by the Kronecker's delta 



(2.10) 



where we have assumed that the a-th bond at site i and 
/3-th bond at site j are to be connected, and that the ten- 
sor at site j is of the modified type at the /3-th position. 



k) = l(M))> for ir = E a V 



(2-3) 
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B. Vertex-Model Interpretation 



and for S = 2 square-lattice VBS model, we have 



We can regard the tensor 4(<t|{?7}) as the Boltz- 
mann weight of a statistical-mechanical model where 
both the spin variables {cr,} ("vertex spins") and the 
bond variables {rf a } are fluctuating variables, which we 
may call spin-verte x m odel. Then the unnormalized 
VBS wavefunction (|]^) is the partition function of the 
spin-vertex model with fixed vertex-spin configuration 
{(7!, (72, . . . , &n}- 



To obtain the quantum-mechanical expectation (l.f ), 
we should calculate the norm (VBS|VBS) which can also 
be regarded as a partition function of a lattice model, 
namely, the vertex model. Assuming that the local ten- 
sors {Ai^ailrf)} are real, we can write the vertex weight 
for the M-valent site as 

Wlfab), . . . , (r? M , fa)] = E M<r\{Ti})M*\{t})> 

<J 

(2.11) 



The sign factor in (2.6) can also be taken into account 
by introducing the modified vertex weights, for example, 



efe,i)e&,0[-.(lU; 



(2.12) 



In this vertex model, the double index (rj a ,£ a ) (a = 
1,2,.. .) plays the role of the bond variable in the or- 
dinary vertex model. Since each component index takes 
two values ±1/2, the model is a 4-state vertex model. 
The quantum-mechanical expectation ([0]) can then be 
regarded as a statistical-mechanical average in the vertex 
model. 



C. Square and Honeycomb Lattices 

We treat the VBS-state-associated vertex model by the 
transfer matrix method. For ease of construction of the 
row-to-row transfer matrix, we restrict the analysis to 
the two cases in the present paper: S — 2 VBS model 
on the square lattice and the S — 3/2 VBS model on 
the honeycomb lattice. In the former, the row-to-row 
transfer matrix can be constructed in the conventional 
way. In the latter, we can map the model on the 
honeycomb lattice to the one on the square lattice J2(]] 
as shown in Fig. ||. Hence either case can be treated as 
a vertex models on the square lattice. 



fig-2 



From the general formula (2.4), we can explicitly write 
down non-zero components of the local tensor. For S — 
3/2 honeycomb-lattice VBS model, we have 



i4(3/2|l/2,l/2, 1/2) =1, 
4(l/2|-l/2,l/2,l/2)=-L 



(2.13) 



4(2|l/2, 1/2, 1/2, 1/2) = 1, 
4(l|-l/2, 1/2, 1/2, 1/2) =J, 



4(0|-l/2,-l/2,l/2,l/2) 



1 



(2.14) 



Other non-zero components are obtained by using the 
permutational symmetry of 4(cr|{^}) in {rj} and the 
property Q2.5Q . The VBS states constructed from these 
tensors are exact ground states of the Hamiltonian where 
the local nearest-neighbor (nn) -type pair Hamiltonian 
hij is the projection operator into the local spin-2S' space, 
p.2 5 j (5 = 3/2 for honeycomb lattice, and S = 2 for 
square lattice): ]lq,pi| 



(2.15) 



Pff |VBS) = 0, for all nn site pair ij. 



D. "Deformed" VBS-type Model on the Honeycomb 
Lattice 

In the original VBS model, the local Hamiltonian P^ s 
is a function of the inner product SfSj, hence, is isotropic. 
p8| , p9| Anisotropic generalization of the VBS models 
and their associate ground-state vector has been known. 
p2| , ^3| Similar generalization has also been known for the 
honeycomb VBS model. |20| The generalized honeycomb 
VBS model contains a "deformation" parameter which 
corresponds to the xxz-type anisotropy of the Hamilto- 
nian. The associated wavefunction is also of the TPFA 
form with a deformed local tensor 4(cr|{?7}) a where a is 
the deformation parameter. Explicit form of 4(cr|{?7}) a 
is simply given by that of the isotropic case ( 2.13j ) with 
substitution 



71 - 1/a - 



(2.16) 



By ^ / vBS(a)( cr i, a 2i ■ ■ ■ , cat), we denote the TPFA wave- 
function made from {A i (ai\{r]' l }) a }, and by |VBS(a)), 
corresponding state vector. Since |VBS(— a)) is essen- 
tially equivalent to |VBS(a)) (via a S^-diagonal unitary 
transformation), we treat a > in the present article. 
(We have confirmed that calculated results for —a are 
identical to those for a.) As has been shown in Ref. 
p(| there actually exists Hamiltonian with nn-type in- 
teractions whose ground state is precisely |VBS(a)) (x 
normalization factor). 



III. CALCULATION METHOD 

In the previous section, we explained how to translate 
physical quantities from the quantum-mechanical frame 
(Eq. (|0) ) into the cl assica l statistical-mechanical one 
(Eq. ( fHll ) and/or Eq. ( fO^ )). 
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In this section, to deal with the classical statistical- 
mechanical frame, we introduce a transfer matrix. So 
now we focus on a way of renormalization of an asymmet- 
ric transfer matrix by using the "classical PWFRG" jl3) . 
Our transfer matrix is composed of 4-state vertex 
weights. We denote this 4-state vertex weight by 
£?(a, 6), which is depicted by Fig. |[ We regard this 
EP(a, b) as a matrix element of E?. 
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A. Magnetization 

Let us explain how to obtain the expectation value of 
any one-point operator Ai at site i (by the PWFRG): 



(A) 



(3.1) 



where ip represents the ground state of the VBS model 
(for instance, \ijj) = |VBS(a))). We represent the numer- 
ator of Eq. fl3.ip by transfer matrices T and Ta, whose 
definitions will be shown later (Eq. (3.3) and Eq. (3.4)): 



(i,\A i \i;)=Tr[T A (T) N - 1 ], 



(3.2) 









M 




j=-(M-l) 







where N is a linear size. T and Ta in this equation 
are defined in the following (for short , we represent 
{• • -P-i,Po,Pi, ■ ■ ■} by {p}). 



(3.3) 



where 2M is another linear size (along this direction, we 
do not impose a periodic boundary condition). And 
is composed of a quantum-mechanical operator A which 
is sa ndw iched between classical weights (Eq. (2.4) and/or 
Eq. (2.9)). Note that E = (id: a local identity oper- 
ator). Using these notations, 



M 

- cn{8 - pyJS s II ^ 

j=-(Af-l) 

Furthermore T is decomposed into 
T = RDL, 



(3.4) 



(3.5) 



where R is composed of column right eigenvectors, and D 
is diagonal matrix whose diagonal elements {Xj} (these 
Ai, A2, . . • are in descending order of the absolute value) 
is the set of eigenvalues of the transfer matrix T and L is 
consist of row left eigenvectors. Here we should keep in 
mind that T is, generally, an asymmetric matrix, upon 
whose left and right eigenvectors L and R we impose the 
normalization condition 



Using Eqs. (Q and (f]| and (§J), 



N-l 



(LT A R) 1 , 1 (X 1 ) If - 1 
(illAlrXAx)"- 1 , 



(3.7) 



where the subscript "1" represents the eigenstate with 
the largest eigenvalue of T in absolute value and we de- 
note its left eigenvector and right eigenvector by (l\ and 
\r) respectively. 

The denominator of Eq. ( |3.l| ) (the "partition func- 
tion") is obtained by 



<*/#) = Tr [T N ] £* Af 



(3.8) 



From Eqs. (fy), (|3/7j) and (|3 
tion 
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Ai ' 



we get one-point func- 



(3.9) 



LR= 1. 



(3.6) 



To get Ta, |r) and Ai in this equation, we use the 
PWFRG method (infinite method). As we stated ear- 
lier, this transfer matrix is asymmetric. So we decide 
how to renormalize the asymmetric transfer matrix us- 
ing the PWFRG. This is because there are several ways 
to treat an asymmetric transfer matrix. In this paper, 
we treat it in the following way which is different from 
our previous treatment. Q In our previous treatment, 
we adopted asymmetric density matrices in the DMRG. 
Unlike this, in this paper, our treatment is equivalent 
to treating symmetric density matrices in the PWFRG 
(namely, we use the singular value decomposition (SVD) 
for left and right eigenvectors of T) . It is easier to per- 
form the renormalization with the SVD than renormal- 
ization with diagonalization of asymmetric matrices, be- 
cause the former is free from difficulties of manipulating 
complex numbers. In compensation for this difficulty, 
we only have to intro duce identity operators as will be 
shown later (see Eq. (3.11)). These identity operators 
is very important for the correct renormalization of an 
asymmetric transfer matrix. 

First of all, we present the outline of our prescription 
in this paper before explaining details. 

Our prescription in this paper is: 



P-l Multiply left eigenvector and right eigenvector ob- 
tained in the previous iteration by transfer matrix 
(and by an identity operator) several times to get 
the improved left eigenvector and right eigenvector. 

P-2 Perform SVD for the improved left and right eigen- 
vectors, respectively to get "projectors" . 

P-3 Select important states and renormalize transfer 
matrices, identity operators and one-point opera- 
tors. 
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P-4 Get improved "projectors" using the "recursion re- 
lation" and, from these "projectors" , construct 
new left and right eigenvectors. 

P-5 Eval uate one-point functions such as (Ai) by using 
Eq. ( |3.9| ). Check the convergence of both left and 
right wavefunctions and the value of the one-point 
function. If convergence is true then we calculate 
2-point function or evaluate A2 to calculate the cor- 
relation length (see section III B| ) . 



P-6 Return step P-l. 



The followings are the details in the above prescription: 
(Hereafter, we use Greek indices to represent the block 
state.) 

First, a new notation is introduced. We renormalize 
bare T into renormalized transfer matrix T 



a',i',j' ,P' 



E [T L D> (*) [T R ]^, (k), 



(3.10) 



fc=i 



where, in this equation, Tl and Tr are renormalized as 
Eq. ( 3.15| ). Eq. ( [3. 10|) is expressed according to Nishino's 
diagram || in Fig. |J 
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Second, we give caveats in our prescription below: 
As for P-l: 

To concretely explain what is mentioned in P-l, we take 
the case of improving right eigenvector |r)i mp . It is im- 
portant that whenever we multiply by a transfer matrix, 
we simultaneously must multiply by left and right renor- 
malized identity operators idr, and idR. Assuming that 
k is some integer, 



Qp - ((id L ) T (id R ) T T 



(3.11) 



where (id^) is the transpose of idL etc. More con- 
cretely, rewriting this equation in its component (we treat 
the case k = 1 for notation simplicity), 

(M ; M>l r )imp = 

E [id L £ [±dKtT^ !0 ,(a',i',j',l3'\r), (3.12) 

a,0,a f ,i' .j' 

where [idi^]" and [id R ]^ are left and right renormalized 
identity operators (see Eq. ( |3.16 )). 

As for P-2: 
SVD is performed as: 



(Z|a,i,j,/3>=5> L ]^ [<7i] M [V L ] 



3,(3 



(3.13) 



(a > t,i, J 8|r)=X;W Wk], [Vr] 



3,13 



(3.14) 



where and {[<tr] },j are sets of singular values. 

And Ul, Vl, Ur and Vr are orthogonal matrices. 
As for P-3: 

We present examples. New renormalized transfer matrix: 
is given by 



-(new) 



(0= E W W(fc)£5(*,o 

a,f3.k 



(3.15) 

And the left identity operator: idL is renormalized as 



id 



(new) 



13 



(3.16) 



At the very first iteration (when we treat a 4-site chain 
of vertices T = EEEE), we set 



[id L ]fl = S, 



a - °<*,p- 



(3.17) 



This is an initial condition for the recursion relation 
Eq. O 



Next, we prepare an ob ject T Aj,. This is for Eq. (3.21). 
In the same way as Eq. (|3.10| ), we renormalize bare Ta 
into Ta as follows. 



Txjf^p, = E (*) M'S' (*)> ( 3 - 18 ) 

fc=i 

where TAl is made of 

M;!(0= E W i^];(k)[E A ]](k,i) [u R ]f. 

a,/3,k 

(3.19) 

As for P-4: 

For example, Ur, 's recursion relation is: 



(new) 
R 



01,3 



= E W 



(old) 



(3.20) 



As for P-5: 

We show how to obtain the expectation value of any one- 
point operator Ai in the next way. 



(Ai) = 



(l\a,i,j,P)T Aaltild , i01 {a ,i ,j ,/3\r) 

<x,i,j,P,a',i',3',l3' 

(3.21) 



and 
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B. Correlation Length 

In what follows, we explain two methods of calculating 
correlation lengths. For concreteness, we explain the case 
of the spin-3/2 honeycomb lattice model. It is also easy 
to deal with the spin-2 square lattice model. 



1. Calculus by observation of dumping correlators 

We set j — i = r(> 0) and assume N ^> r ^> 1 in the 
sequel. And we locate each operator A in the position 
"b" in Fig. ~ 



So, we have like Eq.(3.7) 



(AiAj) 



= (l\T A [T- l \T A \r){\{)-<r +1 \ 



(3.22) 



From this equation, it follows that we measure corre- 
lation length in the direction depicted in Fig. ||. 



fig-5 



In off-critical region, we expect 
lim (^)-(^)<i4 i )SCexp(-r/0/r w . (3.23) 

j—i — >oo 



We bear in mind that in Eq. (3.22), £ and r includes 
C = 1 and t = 1 (let us recall that the subscript "1" 
represents the eigenstate with the largest eigenvalue of T 
in absolute value). So, in the region (Aj) — 0, we easily 
evaluate correlation length from Eq. ( 3.22 ) and Eq. ( 3.25 ) 
as 



£, 1 = lim log 

j—i — >oc 



(AiA j+1 ) 



(3.24) 



and (Z2I and |ra) are the left and right eigenvectors of the 
eigenvalue A2, respectively 



From Eqs. fl3.25D , ( |3.26| ) and (|3_27j), we evaluate A 2 to 
obtain the correlation length We use the power 

method to evaluate A2 that is, we perform the power 
method in the space which is orthogonal to the ground 
state vectors (l\ and |r), as follows. 

First, we prepare the initial vectors and |fW). 

And then we get the improved left eigenvector (l^ n >\ and 
the improved right eigenvector |fW). which are orthog- 
onal to (Z| and |r). Namely, 



(Z»| = (H n -V\T- 



(|(»-i)|r|r) 



and 



|f< n >) = T\^ n -^) - \r) 



(l\r) 

^| T | f (n-l) 



(3.28) 



(l\r) 

We should note that our (l^] and |fW' 
Eq.(E|). So, slightly modifying Eq. (P 



(3.29) 

do not satisfy 
, we calculate 



2nd largest eigenvalue A2 of T in the following way. 



A, 



lim A 



(») 



where a!™-* is 



(n) _ ([(")|T|f(")) 



Ar J = 



(3.30) 



(3.31) 



It should be noted that we save memory resource in a 
computer by using the power method (That is we need 
only 0(m 2 ) memory, not 0(m 4 ), where m is a number of 
retained base in the PWFRG). So, if necessary, we can 
perform the above calculation with a larger number m. 



IV. RESULTS 

A. Isotropic VBS models on honeycomb and square 
lattices 



2. Evaluation by using the power method 

Another method by which correlation lengths are eval- 
uated is the following. B y ma king use of Eq. (3.5), we 
can further transform Eq.( 3.22 ) into 



(AiAj) S (A^iAj) + C A exp[-r/£ pro ], 

where a constant Ca is 

^ = (l\T A \r 2 ) (l 2 \T A \r) A 2 
A2 A2 Ai 

and the definition of £ pm in this equation is 



(3.25) 



(3.26) 



(3.27) 



As has been mentioned in Sec. |J, both the S = 3/2 
VBS model on the honeycomb and the S = 2 VBS model 
on the square lattice can be mapped to classical vertex 
models on the square lattice. We made PWFRG calcula- 
tion for the sublattice magnetization and the correlation 



length in a way as described in Sec. Ill 



We have confirmed that sublattice magnetization is 
zero (Neel order is absent) in both models; both mod- 
els are in the disordered phase, in agreement with previ- 
ous studies. Jl||,|l9],|4| As for the correlation length £, we 
obtain 



r 1 



honeycomb 

r 1 

S square 



1.67 
0.52 



(4.1) 

We should note that, in the honeycomb case, the length 
unit for the above value of £hon eyC omb ^ s ^ e ^ a ^i ce spacing 
of the mapped square lattice. 
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B. S = 3/2 deformed VBS model on the honeycomb 
lattice 



Unlike the isotropic case, the anisotropic (xxz-lihe) 
generalization of honeycomb VBS model can be either 
in the disordered phase or the ordered (Neel) phase, de- 
pending on the "deformation parameter" a. pQ ] 

The PWFRG allows us to study the behavior of the 
system in much more detail, which we present in this 
subsection. In the actual calculation, we have observed 
that the PWFRG calculations are rapidly convergent in 
increasing m (number of retained bases), except for the 
cases very near the transition point. Consequently, for 
most values of a, very small to, say m = 12, is sufficient 
for the accuracy required in the present study. 



1. Staggered magnetization 

In Fig. ||, we show the PWFRG result of the sublat- 
tice magnetization (staggered magnetization) M st as a 
function of the parameter a. 



fig-6 



The Neel order which exists in the large-a region dis- 
appear at a — a c = 2.54, in agreement with the Monte- 
Carlo result. J2(| It may well be expected that the Ising- 
like anisotropy make the system be in the Ising universal- 
ity class as regards this anisotropy-induced phase transi- 
tion. To check this, we plot (M st ) s versus a as shown in 
Fig.j^. The clear linear behavior near a c implies 



ng-7 



M st ~ {a-a c fi 
with £ = 1/8. 



2. Correlation length 



(4.2) 



In the transfer-matrix method, the correlation length 
is usually calculated from the ratio between the largest 
and the next-to-largest eigenvalues of the transfer ma- 
trix. This definition of the c orrelation length is, in our 



shown in Fig.||. 



notation, £ pm (see Eq. ( p.27|) ), whose PWFRG result is 



fig-8 



The linear behavior of £T]^ near a c is consistent with 
the expectation that the system is in the Ising-model 
universality-class, namely, 



r 1 

Spm 



with v 



(4.3) 



In addition to the critical behavior at a c , a notable fea- 
ture of £ pm is that it exhibits cusp singularity at the 
isotropic point a — a[ BO = \/3. To clarify the origin of 



this cusp, we calculated longitudinal correlation length 
£i(a) and the transverse correlation length £t(a) from the 
correlation functions (SfSj) and (SfSJ), as shown in 
Fig|_ 



fig-9 



Totally different a-dependence between ^ 1 (a) and 
£ t -1 (a) is seen; the isotropic point a = a/3 is the crossing 
point of the two curves. Since £ pm is determined by the 
correlator with the largest correlation length (or, equiv- 
alently, smallest inverse correlation length), we have 



f -i = J Ci _1 (a) (o> V3), 
v "" 1 £t -1 (a) (o<V3), 



(4.4) 



which explains the cusp at a = a- lso — a/3 in Figj^. 

It is interesting to compare the correlation lengths in 
the 2D deformed VBS model with those in the ID coun- 
terpart (deformed 5 = 3/2 chain ) |p3| , which we denote 
^ (1D) (a) and £ t (1D) ( a )- Explicitly, we have Q 



= log 



V(a 2 + 


l) 2 + 8 


| a 2 


-1| 


V(a 2 + 


l) 2 + 8 



(4.5) 



(4.6) 



which also show the crossover behavior at the isotropic 
point a = a/3- It should be noted that magnitudes of 
these ID correlation lengths are almost the same as those 
in the 2D case in a wide range of parameter values. Major 
differences are the absence of the critical point a c and the 

ID) 



non-vanishing behavior of f^ 11 
(s ee Fig. Ii0| ). 



fig-10 



the 



limit 



In Ref . |20(] , the "asymptotic equivalence" between the 
deformed honeycomb VBS model in the large-a region 
and the free-fermion model is pointed out, which gives 
fairly accurate estimation of the critical point a c . As a 
further check of this equivalence, we compare the cor- 
relation length £i with that of the free-fermion model 
which we denote by £ff = £ff(a). The method given in 
Refs. [p5|-p7| allows us to obtain £ff(a) from the solution 
of Q(iuj*) = 0, where Q(<j>) is given by Eq.(24) in Ref. |2r} 



in the disordered phase and £ ff 1 = 2u>* in the 



ordered phase). In Fig.[TT], we see a remarkable agree- 
ment in a unexpectedly wide parameter range including 
the small-a region where the "asymptotic equivalence" 
may not hold any longer. 



fig-11 



The cusp behavior of the correlation length implies 
that the isotropic VBS model may be a "singular" point 
in the parameter space. We should point out that a sim- 
ilar cusp-like behavior of the correlation length at the 
VBS point has been known in the /3 — £ -1 curve j2§] and 
the 9 — £ curve fl29] of the S = 1 bilinear-biquadratic 
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chain, where —(3 and 9 relate to the relative amplitude of 
the biquadratic exchange term. To explore the relation 
between these cusps in the different parameter spaces is 
an interesting problem, which may help us to character- 
ize the disordered phase in the 2D VBS-like models, in 
more detail: like the string order parameter which 
characterizes the Haldane phase in the ID case. 



V. SUMMARY AND CONCLUSION 

In this paper, we have presented a new way of applica- 
tion of the product wavefunction renormalization group 
(PWFRG) which is a variant of the density-matrix renor- 
malization group (DMRG), to two-dimensional (2D) 
valence-bond-solid (VBS) type models. The exact 
tensors-product-form (TPF) structure of the ground- 
state wavefunctions have allowed us to map the sys- 
tems into two-dimensional classical statistical-mechanical 
models, namely, the vertex models which can be treated 
by the transfer-matrix method. We extend the PWFRG 
method so that it can handle the asymmetric transfer 
matrix associated with the VBS-type models. 

For the isotropic 2D VBS models on the honeycomb 
lattice (S — 3/2) and the square lattice (S = 2), we have 
confirmed that they are in the disordered phase without 
the Neel order. We have obtained accurate values of cor- 
relation length for both models. For the anisotropically 
generalized (or deformed) VBS-type model on the honey- 
comb lattice pq| , we have obtained detailed parameter- 
dependence of the Neel order and the correlation length. 
We have confirmed the anisotropy induced phase tran- 
sition at a critical value of the deformation parameter, 
and that this second-order phase transition belongs to 
the 2D Ising-model universality class. In the disordered 
phase, we have found a cusp-like behavior of the corre- 
lation length. We have explained this behavior as the 
crossover between the longitudinal and the transverse 
correlation length. Further, the PWFRG-calculated lon- 
gitudinal correlation length well agrees with the corre- 
lation length of the free-fermion model in a fairly wide 
range of the deformation parameter, which far exceeds 
the prediction of the "asymptotic equivalence" . [|2(J 

As has been mentioned in the Introduction, the present 
study is also aimed at higher-dimensional generaliza- 
tion of the DMRG-type numerical renormalization ap- 
proach, based on the TPF-ansatz on the trial wavefunc- 
tions. From this point of view, what we have done in the 
present paper is calculation of physical quantities under 
given TPF-wavefunction with fixed local tensor. What 
we should do next is to vary the local tensor to make 
the variational calculation. A "direct numerical" varia- 
tional calculation can actually be done, where we sweep 
the variational parameter and numerically find the min- 
imum of the energy-expectation. The PWFRG method 
is also useful for reliable determination of the optimal 
variational value of the parameter. Studies along this 



line is now undertaken, whose details will be published 
elsewhere. 
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FIG. 10. Comparison of for 1-dimensional (S = 3/2) 
case and one for 2-dimensional case (S = 3/2) in the dis- 
ordered region. Circles (o) and squares (□) are represented 
£,^ 1 (a) and ^ 1 (a), respectively. A solid curve is (£[ 1D \a))~ 1 
and broken curves are (^ 1D ' > (a)) -1 • 



FIG. 11. Inverse correlation lengths of the free-fermion 
model (solid curves) and honeycomb VBS model (o:£j~ (a), 
+:£pm(a))- 



FIG. 1. The pictorial representation of A((j|{?7}) of a 
M-valent site. 



FIG. 2. The mapping of the model on the honeycomb 
lattice to the model on the square lattice. Circles (•), which 
express quantum spins (S = 3/2), are connected by solid 
lines that form the honeycomb lattice. The small letters "a" 
and "b" ind icate sites of the two sublattices, respectively (see 
Sec. pIBl\). 



FIG. 3. Eq(a,b) is expressed according to Nishino's dia- 
gram [H] . Each circle (o) represents the mapped classical site, 
whose states are expressed by "p", "q", "a" and "b". These 
state variables take 4 states of the double index (jjJ,,^). 



_FIG. 4. 

± a' ,i> ,j',0'- 



The renormalized transfer matrix element 



FIG. 5. The direction along which correlation lengths of 
the model on the honeycomb lattice are measured (a solid 
line) . 



FIG. 6. The PWFRG result of the sublattice magnetiza- 
tion (staggered magnetization) M st as a function of the pa- 
rameter a. 



FIG. 7. (Mst) 8 
a c (— 2.54) is seen. 



versus a. The linear behavior near 



FIG. 8. as a function of a. 



FIG. 9. The inverse of the longitudinal correlation length 
^j - (o) and the inverse of the transverse correlation length 
^ t ~ 1 (□), in the disordered region. 
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